Colliding clouds of strongly interacting spin-polarized fermions 
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Motivated by a recent experiment at MIT, we consider the collision of two clouds of spin-polarized 
atomic Fermi gases close to a Feshbach resonance. We explain why two dilute gas clouds, with 
underlying attractive interactions between its constituents, bounce off each other in the strongly 
interacting regime. Our hydrodynamic analysis, in excellent agreement with experiment, gives 
strong evidence for a novel metastable many-body state with effective repulsive interactions. 
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I. INTRODUCTION 

Ultracold atomic gases [l], open up new frontiers in 
the study of non-equilibrium dynamics of strongly inter- 
acting quantum systems. In contrast to electronic sys- 
tems, the Fermi energy in quantum gases is of the order 
of a few kilohertz, so that questions about metastability 
and equilibration in quantum systems can now be ex- 
plored in real time in the laboratory. Exploring these 
new regimes and developing new theoretical tools to gain 
insight into the non-equilibrium dynamics of many-body 
systems is an exciting challenge at the intersection of 
condensed matter and atomic-molecular-optical physics. 

The recent experiment of Ref. Q offers a beauti- 
ful example of non-equilibrium dynamics that is com- 
pletely unexpected and, at first sight, counterintuitive. 
Two clouds of ultracold fermions prepared in different 
hyperfine-Zeeman ("spin") states a — f and 1 are ini- 
tially separated using a Stern-Gerlach field; see Fig. Q] 
Once the field is turned off, the clouds are impelled to- 
gether by the confining harmonic potential and they col- 
lide. In the strongly interacting unitary regime these 
very dilute clouds are observed to bounce off each other, 
almost as if they were colliding billiard balls! This is 
truly remarkable in view of the fact that the underly- 
ing atomic interaction is attractive, and the equilibrium 
ground state is known to be a paired superfluid. Why 
then do the clouds behave as though the interactions are 
repulsive? 

Two additional aspects of Ref. [3j are also noteworthy. 
Once the short-time bounce is damped out, the centers 
of mass of the clouds remain separated for a considerable 
duration at intermediate times (~ 100 ms), before even- 
tually merging. In fact, the analysis in Ref. 0] focused on 
precisely the slow merging of the clouds, or spin diffusion, 
at long times (~ 0.5 s). Finally, the observed dynamics 
changes qualitatively as one moves out of the strongly 
interacting regime, with the clouds merging rapidly for 
weak interactions. 

In this paper we gain insight into three aspects of this 
remarkable dynamics using a hydrodynamic approach in 
the strongly interacting regime. (1) We explain why the 
clouds repeatedly bounce off each other at short times. 
(2) We also explain why the centers of mass of the two 



spin up cloud 



spin down cloud 




Harmonic Trap 

FIG. 1: Schematic of the colliding cloud experiment. 

clouds remain separated at intermediate time scales. We 
show that these results are natural consequences of a 
metastable many-body state on the "upper branch" of 
the Feshbach resonance where the effective interactions 
are repulsive 0, (3) Finally, we provide insight into 
how the dynamics changes as a function of the strength 
of the interaction. 



II. HYDRODYNAMICS 

Our analysis is based on two hypotheses: (i) The col- 
lisions between atoms are sufficiently rapid to establish 
local thermodynamic equilibrium when the clouds over- 
lap, (ii) The loss from the upper branch scattering state 
is slow in the experiment of Ref. 3]. At unitarity, the 
two-body collision rate is very large, l/r% ~ €f (with 
h = Ub = 1), for a range of temperatures 0.1 < T/ep < 
0.3 @, 0]- With a Fermi energy ep ~ 10 4 Hz and an 
axial trap frequency lo z ~ 10Hz H, the gas in the over- 
lap region of the clouds will reach local thermodynamic 
equilibrium within ~ 10~ 3 trap periods. Thus the clouds 
behave hydro dynamically in the overlap region Q and 
their dynamics reflect the (metastable) equation of state 
of the gas. (Hydrodynamics also describes well the very 
different behavior observed in colliding clouds of spin- 
balanced superfluid gases at unitarity [9j.) 

As the two clouds come together in the overlap re- 
gion, the atoms are initially in scattering states, and the 
formation of two-body spin-singlet bound states requires 
three-body collisions in order to satisfy kinematic con- 
straints. The time scale r 3 for such processes is not well 
understood for strongly interacting gases, however, there 
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are indications that it is enhanced close to unitarity ■ 
In addition, the r 3 relevant to Ref. Q is further enhanced 
relative to the microscopic decay time since three-body 
collisions arc limited to a small overlap region for a frac- 
tion of the oscillation period 

Now, for the experimentally relevant time interval 
T2 <C t <C T3, the rapid two-body collisions have already 
established thermal equilibrium, while the slow three- 
body recombination process has yet to drive the system 
into the lower branch. This allows us to model the cloud 
dynamics using Euler's equations 
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with an appropriate energy density f [^(r)] for the scat- 
tering states (see below). Here v CT and n a are the velocity 
and density of the a fermions, Vt rap = m{Lo\p 2 +w 2 z 2 ) /2 
is the trap potential. We include a phenomenological 7 to 
account for st rong sp in-current damping in the hydrody- 
namic regime [111 Il2j . Viscosity, describing the damping 
of in-phase current (v-j- = V|,), is ignored since the motion 
of the colliding clouds is primarily out-of-phase @, [HI . 

As a result of tight confinement in the radial direc- 
tion, the two modes with the lowest energies are the 
spin-dipole mode and axial breathing mode. Hence, we 
concentrate on the centers-of-mass of the clouds and their 
widths along the z-direction and reformulate (HJ and @ 
in terms of those variables. Using the notation (■ ■ • ) a = 
1 /Na J d 3 rn a (- • ■ ), we define the center-of-mass positions 
z a = (z) a and widths Sz a = \/&((z — z a ) 2 ) a . The \f% 
ensures that bz a coincides with the axial Thomas-Fermi 
(TF) radius R z in equilibrium. ((T|) and © then lead to 
the exact equations of motion 
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The problem can be further simplified while still re- 
taining the essential physics of the (nonlinear) coupling 
between the spin-dipole and axial breathing modes by 
using a TF ansatz: 
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where R a = y / 2e t p/muj 2 is the TF radius along the a- 
axis and e° F = (ui 2 1 _uj z ) 1 /' i (2>N) 1 /' i is the chemical poten- 
tial of an ideal two-component Fermi gas (iV-j- = N\, — 
N/2). This ansatz allows for a time-dependent center- 
of-mass Z(j as well as axial compression of the clouds. 
The continuity equation ((2J leads to the velocity field 

Vcr = V a Z With V a (z,t) = Z a — ZcrSZa/SZv + zSz a /Sz a . 

Axial and radial symmetry lets us set z\ = z = — zj, 
and 5z^ — 8z± = Sz. Using <j5j in ((3]) and (j4|) gives the 
coupled nonlinear integro-differential equations: 
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In the unitarity regime, the spin-current damping 12 1 
assumes a particularly simple form when the relative ve- 
locity is of order of the Fermi velocity, as is the case at 
early times for two clouds initially separated by R z . In 
the overlap region, we take 



7 = Jy/ e F-t e Fl' 



(8) 



where 7 is of order unity at unitarity and eir CT (r) = 
(67r 2 n cr (r, i)) 2 / 3 /2m is the local Fermi energy. This sim- 
ple choice of damping is, if anything, an overestimate, 
and we have checked that our results are robust against 
reasonable modification of the damping parameter. 



III. ENERGY FUNCTIONAL 

We now need to specify the energy functional £[n a ] 
relevant for dynamics at times Ti <C t <C 13. The ground 
state ( "lower branch" of the Feshbach resonance) for scat- 
tering length a s > necessarily involves bound pairs. 
But for t <C T3, the three-body processes required to re- 
lax to this state have not yet occurred. The system has, 
however, developed two-body correlations characteristic 
of the metastable "upper branch" state, studied theoret- 
ically in Refs. 0, H, EH, motivated by an earlier experi- 
ment [l5j • The approximate many-body wavefunction in 
the upper branch is 
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where $s({ri cr })'s are Slater determinants and the Jas- 
trow factor f(r) describes the short-range correlations 
between fermions. The effective repulsion between ^ and 
1 atoms in the upper branch is crucially related to the 
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FIG. 2: (Color online) Time dependence of the center-of-mass 
positions z CT of the two spin clouds at unitarity, with the red 
(upper) and blue (lower) curves denoting the two spin species. 
We use the upper branch energy functional (see text), with 
corresponding to N = 7.5 x 10 5 atoms, trap anisotropy 
u)±/ujz — 10, and damping 7 = 1. Initially the two clouds are 
displaced from the trap center by the axial Thomas-Fermi 
radius, %(t = 0) = R z , %i(t = 0) = — R z with initial axial 
width Sz(0) = R z equal to its equilibrium value. 



node in f(r), in contrast to the nodeless Jastrow fac- 
tor for the lower branch. For small a s > 0, the node 
occurs at a s , with f(r) ~ (1 — a s /r) similar to the two- 
body problem. For large a s , the node saturates [j| to 
~ 1/kp, the only length scale at unitarity. Quantum 
Monte Carlo (QMC) studies 0, U of the upper-branch 
wavefunction (|9]) reveal that the system undergoes (fer- 
romagnetic) phase separation for sufficiently strong in- 
teractions kFCLs ^ 1- 

We need £ [n CT ] for arbitrary polarization, which has 
not been studied by QMC. We thus use the lowest or- 
der constraint variational (LOCV) approximation [l6j . 
which has been used for the upper branch [l7[ and is in 
close agreement with QMC data [3, H[ . As shown in the 
Appendix, the upper branch LOCV energy density is 



£ = -ejrfnt + -epini + ttOt^ + %^t)- ( 10 ) 



The interaction energies A^ and A^ are functions of fcj?a s 
and x = ni/n-f. Hydrodynamics requires T > 0.1 ep Q; 
for simplicity, we use the zero temperature LOCV energy 
functional to study the dynamics, expecting it to quali- 
tatively describe the physics at low temperatures. It is 
important to note here that while we rely on the approx- 
imate LOCV calculation to obtain the equation of state, 
the underlying physics of the effective repulsive interac- 
tion is independent of the particular scheme used and a 
more precise calculation would only leads to quantita- 
tive improvement without changing the crucial physical 
picture of the bounce. 
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FIG. 3: (Color online) Time dependence of the center-of-mass 
positions z a of the two spin clouds for fcF(0)a s = 2. The 
parameters are the same as in Fig[2] 

IV. DYNAMICS AT UNITARITY 

We solve © and Q, using dTUJ) and ©; for details 
see Appendix [Bj The results at unitarity are shown in 
Fig. [2] The two clouds bounce off each other for sev- 
eral oscillations due to the repulsive nature of the upper 
branch functional. The period of the initial bounce is 
roughly 0.56(27r/w 2 ), slightly less than the experimen- 
tal value 0.61(27r/o; 2 ) 0. We attribute this difference to 
the simple ansatz ([5]) used to model the dynamics. The 
bounce persists for several cycles in spite of the very large 
damping (c f /uj z ^> 1) because the overlap between the 
two clouds is small. Once the oscillation is damped out, 
the centers of mass of the two clouds remain segregated, 
with a final {z^ — z^) ~ 0.4 times the initial separation. 
This intermediate time behavior reflects the tendency for 
system to phase segregate in the upper branch at unitar- 
ity. 



V. COLLISIONAL DYNAMICS AWAY FROM 
UNITARITY 

Away from unitarity, it is harder to reach the hydrody- 
namic regime. If u z is sufficiently small, however, it will 
always be the case that the dynamics in this direction 
are hydrodynamic. The greater challenge at finite a s > 
is that T3 may not be much greater than t^- As a s de- 
creases from unitarity, one expects that T3 reaches a min- 
imum for kpa s ~ 1 and then becomes large again [lg| . 
with T3 ~ (nal)~ 2 ep 1 for kpa s -C 1. In contrast, the 
cross-section a 2 s determines T2 ~ (naD^^tp 1 . Hence, 
for small kpa s , again T3 3> T2. However, when kpa s ~ 1, 
it is conceivable that T3 ~ t-i- 

With these caveats in mind, we solve (see Appendix [B] 
for details) our hydrodynamic equations away from 
unitarity to understand the relationship between the 
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FIG. 4: (Color online) Time dependence of the center-of-mass 
positions z a of the two spin clouds for fcF(0)a s = 0.5. The 
parameters are the same as in Fig[2] 

intermediate-time dynamics, after the bounce is damped 
out, and the upper branch equation of state: Fermi liq- 
uid (fc,Fa s < 1) versus phase separated (fcfO s > 1). The 
solutions of ([6]) and (O at kp(0)a s = 2 and 0.5 (using the 
same initial conditions and e F as in Fig. [5]) , are shown 
in Figs. [3] and |U Here kp(0) = y/2me° F is the ideal gas 
Fermi wavevector at the cloud center. Appropriate for 
the smaller interaction strength, we use an expression 
for the spin-current damping obtained from kinetic the- 
ory [ill ]. When symmetrized in the spin components, it 
takes the simple form |19| 

7 = (4/97r)( A /A; Ft fc F4 .a fl )VeFten- ( n ) 

Even at small kp(0)a s , the clouds exhibit a weak 
bounce due to compressional recoil, which damps out 
very quickly. We see, however, a clear difference between 
Fig. [3] (kpa s > 1), where the centers of mass remain sep- 
arated, and Fig. \Q(kpa s < 1), where they merge. This 
behavior is qualitatively similar to the experimental re- 
sults shown in Supplementary Fig. 1 of Ref. [3J . 



The long time behavior, not described here, will be dom- 
inated by spin diffusion Q at non-zero temperatures and 
three-body processes. 

Finally, we comment on the experiment of Jo et al. [H[ 
in which a spin-balanced mixture, initially at a small 
a s > 0, is swept close to resonance, generating strong 
interactions in the upper branch. It appears, however, 
that rapid losses render it unstable to the lower branch 
within a very short time (r3 im i cro ~ 1ms) [l(l[2(|[2l| and 
phase separated ferromagnetic domains are not observed. 

In contrast, the specific initial configuration in the col- 
liding cloud experiment Q proves to be crucial for the 
metastability of the upper branch. Three-body loss is 
limited spatially within the overlap region between dif- 
ferent spins at the trap center, and temporally to a frac- 
tion of the oscillation period. As such, the effective t% 
is greatly enhanced (ts 3> T3 tm i CIO ) and one can study 
the metastable up per branch. Note that while T3 im j cro 
is small, ~ 10/ ep [21j, it is still an order of magnitude 
greater than ti ~ ep, and thus local thermodynamic 
equilibrium can be established in the upper branch. 

In summary, we have investigated the short-time 
bounce dynamics of the recent MIT experiment Q. The 
bounce frequency we found is in good agreement with 
the experiments. Furthermore, the experimental obser- 
vation that the two clouds sit side-by-side for as long 
as ~ 100ms validates our assumption of a long T3. We 
argue that the intermediate time behavior of the two col- 
liding clouds indicates that the system favors a phase 
separated metastable state for larger value of kpa s . We 
cannot see how a lower-branch energy functional, with 
attractive interactions forming bound pairs, could lead 
to the observed dynamics. 

Note added: Recently, a Boltzmann approach was used 
to investigate the same problem in the high-T regime [22j . 
The system is found to be strongly hydrodynamic for 
kpa s > 1. 
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VI. DISCUSSION 

We now compare our results in detail with experi- 
ments. As seen from Fig. [2j our results at unitarity - 
the short-time bounce, the damping of the oscillations 
and the separation of the clouds at intermediate times 
- are all in very good agreement with Ref. [3j. Away 
from unitarity, at finite values of kpa s > 0, our re- 
sults are qualitatively similar to the data in Supplemen- 
tary Fig. 1 of Ref. [3J. For times uj^ 1 <C t <C T3, say, 
t = 200 ms, the system remains phase segregated for 
large kpa s [Figs. l(f,g)] s while completely mixed for small 
fcj?a s [Figs. l(c,d,e)]. In fact, a very rough estimate for 
the critical kpa s can be read off from the experimental 
data: it is between 0.26 [Fig. 1(e)] and 1.2 [Fig. 1(f)]. 
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Appendix A: Upper branch and LOCV 

In this Appendix, we present a detailed discussion of 
the so-called "upper branch" of a Feshbach scattering res- 
onance in two-component Fermi gases 0, H[ as well as the 
lowest order constraint variational (LOCV) method [l6| . 
which we use to calculate the upper branch equation of 
the state [13] • 
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1. Upper Branch 

The concept of the upper branch of a Feshbach res- 
onance comes from the two-body problem (either in a 
finite box or in a harmonic potential), where it is com- 
pletely well-defined for all values of a s , positive or nega- 
tive. The two-body wavefunction in the upper branch is 
a scattering state with a single node that makes it orthog- 
onal to the ground state (lower branch) wave-function. 

This notion has been generalized to the many-body 
case [H, H| by writing a Jastrow- Slater wavefunction, 
where the Jastrow correlation factor has a node; see 
Eq. (9) of the main paper and the discussion immedi- 
ately following this equation. 

To summarize our definition of the upper branch in the 
many-body problem, any sensible definition must at least 
satisfy the following conditions ||: 

(1) The many-body wavefunction includes, apart from 
the nodes introduced by the Pauli principle, one addi- 
tional node for any pair of fermions with opposite spin. 

(2) The wavefunction should reduce, in the limit of the 
two-body problem, to that of the scattering states with 
one node in the relative wavefunction. 

(3) The energy of the system must be larger than that 
of the non-interacting Fermi gas and, it should reduce to 
the perturbative result in the weakly interacting regime 
< k F a s < 1. 

We emphasize that orthogonality with the many-body 
ground state (of the BCS-BEC crossover) is not sufficient 
to be on the upper branch. 



Variation of the energy (IA1[) with respect to f(r), while 
taking into account the constraint (|A2|) by a Lagrange 
multiplier A, gives us the "Schrodinger" equation 



v(r) 



f(r) = A/(r). 



(A3) 



Retaining only the s-wave part of this equation, we find 
the general solution with one node is given by 



al sin(/«(r — b)) 
r sm.(n(d — b)) 



(A4) 



where k = \/mX. We find that f(d) — 1 and f'(d) = 
lead to nd = tan(/t(d — b)) and the Bethe-Peierls bound- 
ary condition gives na s — tan nb. With the normalization 
condition (|A2|) . we can solve for k, d and b for each value 
of scattering length a s . The energy of the system is given 
simply by 



1 



-nX. 



(A5) 



In the case of our interest, however, the system is not 
necessarily balanced, so that x = ri\/ri+ need not be 
unity. We have to extend the LOCV calculation to the 
spin-imbalanced case in which the Slater determinants 
will have different sizes. We find the energy of the sys- 
tem is given by 



2. LOCV 

LOCV is an approximation scheme to evaluate the en- 
ergy of the Jastrow-Slater state that originated in nu- 
clear many-body physics and has been used for var- 
ious quantum fluids. Within LOCV, the energy of the 
Jastrow-Slater state [Eq. (9) of paper] is given by 



£ = le F n^^Jd\r(r) 



1- vyr) 

m 



f(r), (Al) 



where v(r) is the short-range two-body potential, n — 
Tif + n_i is the total density, and ep = (37r 2 n)/2m is the 
Fermi energy. The effect of a zero-range contact potential 
may be written in terms of the Bcthc-Pcierls boundary 
condition lim r _>.o ( r f( r ))' I (, r f( r )) — ~ l/ fl s- In addition, 
within LOCV, the Jastrow function f(r) satisfies the fol- 
lowing conditions: f(r > d) = 1 and f'(d) — 0. The 
"healing length" d in turn is defined so that 



drr 2 f 2 (r) = 1. 
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Now, in general, there is no unique way to enforce the 
normalization conditions as in (|A2j) . A natural extension 
is to use both normalizations 



47r(n-?v) / drr 2 fl{r) = 1, 



(A7) 



which introduces two healing lengths, d\ and d^. Accord- 
ingly, we shall introduce two Lagrange multipliers and 
A^ and the energy of the system can then be written as 
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(n t X i + n i X t ). (A8) 



6 



3.5 
3.0 
I 2 - 5 

0) 

2.0 
1.5 
1.0 



^1/fc F a s =0 


(a) 


~-^Wt F a s =2.0 




1//< F a s =10____ 





0.4 0.6 
P 



4.0 
3.5 
3.0 
' 2.5 
2.0 
1.5 
1.Q 







1.0 


* — . — . — . — . — . — . — — 






\ P=0.1 




0.8 
0.6 




\ 

\ 


(c) 






0.4 
0.2 


FP 


PP 




P=1.0 




0.0 









1.0 

Vk F a s 



1.0 

Vk F a s 



1.5 



2.0 



FIG. 5: The single particle energy in the upper branch, obtained from lowest order constraint variational method. The energy 
are in unit of 3/5eF (a), energy as a function of the polarization P for various values of kpa s as indicated in the figure, (b), 
energy as a function of the interaction parameter fc^ffls for various values of polarization P as indicated in the figure, (c), phase 
diagram of the interacting two-component Fermi gases. The two vertical red dashed lines indicates the separation between 
three phases. For l/(kFa 3 ) > 1.124, the system is in the Fermi liquid (FL) phase. For 1.124 > l/(kpa a ) > 0.91, the system is 
in the partially polarized phase (PP) and for l/(kpa s ) < 0.91, the system is fully polarized. 



The above energy functional reduces to the usual per- 
turbative result in the weak coupling limit and also, as we 
shall show later, reproduces the phase diagram of the in- 
teracting fermion system given by Monte Carlo methods 

0, [E[ with very good agreement. We note that the same 
energy functional (IA8|) was analyzed recently by Heisel- 
berg jT3|, in the context of the experiment of ref. [l5j . 

In Fig. [5] (a), we show the upper branch energy (|A8I) as 
a function of polarization P = (1 — x)/(l + x) for various 
values of kpa s . For small kpa s (in fact, kpa s < 0.91), 
the minimum energy is attained at P = 0. For larger 
values, kpa s > 1.124, the minimum value is attained at 
P = 1. For intermediate values of kpa s , the minimum 
value occurs at a value of < P c < 1. In Fig. [5] (b), we 
show the upper branch energy as a function of 1 /kpa s for 
various values of the polarization P. Note that for P = 

1, i.e., completely polarized, the energy is completely 
flat, since there is no interaction between the polarized 
fermions. 

In Fig. [5] (c), we show the phase diagram of the sys- 
tem as obtained from LOCV. For l/(kpa s ) > 1.124, the 
system is a homogeneous mixture of the two hypcrfinc- 
Zeeman states. We call this a Fermi liquid (FL) state. 



For 1.124 > l/(kp a s ) > 0.91, the system is in the par- 
tially polarized phase (PP). Here the transition as pre- 
dicted by the LOCV method is second order and accom- 
panied by a divergent spin susceptibility [13] ■ Note that 
the value of l/(kpa s ) at the transition is very close to 
that predicted by Monte Carlo calculations (3, [El and 
compares favorably with the calculation in Ref. [23| . 
Lastly, for l/(kpa s ) < 0.91, the system is a fully po- 
larized, non-interacting Fermi gas. 



Appendix B: Numerical details 

In this section, we explain how we use the LOCV en- 
ergy in the numerical solution of the hydrodynamic equa- 
tions (6) and (7) of the main text. The important quan- 
tity that enters both equations is 

^ = ^ + 2^ + n ^ +n ^)- (B1) 
We express this in dimensionless form as 



1 d£ 




dX i MA, 2 ~ 1 d\A 

: ^ - + - ) 



r 



(B2) 



Here £ = l/(fc_pa s ) and all energies are scaled by cf^, 
so that A^ = X^/ep^ etc. We note that at unitarity, 
£ = 0, and the expression involves only a single deriva- 
tive dX^/dx that has to be evaluated numerically. For 
arbitrary a s , however, the expression involves two addi- 



tional numerical derivatives. To simplify the numerics at 
"small" a s , we find that we can make a polynomial fit 
to the equation (|B2|) as a function of x and we find 
that, for the values of kp^a s < 2 we are considering, a 
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good fit is provided by the expression where the coefficient 0.27 is very close to the value ob- 



tained by Galitskii [2J, ^(11 - 2 In 2) » 0.259. 

ju t = 1 + ar ( J^t a s + 0.27(fc Ft a s ) 2 J + • • • (B3) 
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